{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.io import mmread\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import scipy.sparse as sp\n",
    "\n",
    "import scipy.sparse.linalg as spla\n",
    "\n",
    "K = mmread('fidap037.mtx') # COO format\n",
    "\n",
    "K = K.tocsc() # for factorization later\n",
    "\n",
    "nR,nC = K.shape\n",
    "\n",
    "b = np.ones((nR,2))\n",
    "\n",
    "b[:,1] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3565, 3565)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "67591"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K.nnz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARAAAAD8CAYAAAC/+/tYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXt4VdWZ8H8vSbijhAiKCJPAoEbFIkLgAeuVghC8dKZeS5vBfgOflU6ZKVMj2oYU0VSxZWpnGPRTSsdSRWY6jcao4KW0UkCwKUKjcknKVYJAlDtJWN8fZ++wz8ne++xzSXKSvL/nOU/OWXvttdYmOS9rvVcxxqAoihIPnVp7AYqitF1UgCiKEjcqQBRFiRsVIIqixI0KEEVR4kYFiKIocZOyAkREbhaRj0Vkm4gUNsP41SLyoYhUiMgGq62PiKwUka3Wz0yrXUTkZ9ZaNonIiBjnel5EakRks6Mt5rlEpMDqv1VECuKcd66I7LGeu0JEJjuuPWTN+7GITHS0x/y7EJGBIvKOiFSKyBYR+W5LPLfPvM3+3CLSVUTWi8ifrbmLrfYcEVlnrf8lEelstXexPm+zrmdHW1OM8/5CRKoczzw8mf/WABhjUu4FpAHbgcFAZ+DPwGVJnqMaOC+i7Qmg0HpfCPzYej8ZKAcEGAOsi3Gua4ERwOZ45wL6ADusn5nW+8w45p0LzHbpe5n179wFyLH+/dPi/V0A/YER1vtewCfWHM363D7zNvtzW2vvab3PANZZz7IcuNtq/0/gfuv9t4H/tN7fDbzkt6Y45v0F8DWX/kn7G0vVHUgesM0Ys8MYcxp4EbitBea9DVhqvV8K3O5o/6UJsRboLSL9gw5qjFkNHEpwronASmPMIWPMYWAlcHMc83pxG/CiMeaUMaYK2Ebo9xDX78IYs88Y84H1/ghQCQxo7uf2mbfZn9ta+1HrY4b1MsCNwAqPZ7b/LVYAN4mI+Kwp1nn9njkpf2OpKkAGALscn3fj/0cQDwZ4U0Q2ish0q+18Y8w+CP0hAv2acT2xzpXMNcy0tq7P20eI5pzX2ppfReh/xhZ77oh5oQWeW0TSRKQCqCH0BdwO1Bpj6l3GaZzDuv45kBXP3JHzGmPsZ55vPfNPRaRLsp85VQWIuLQl2+d+nDFmBDAJeEBErm3l9USbK1lrWAQMAYYD+4CnmnNeEekJ/DcwyxjzhV/XZM7vMm+LPLcxpsEYMxy4iNCuIddnnKTNHTmviFwBPARcCowidCx5MNnzpqoA2Q0MdHy+CNibzAmMMXutnzXAbwj9svfbRxPrZ00zrifWuZKyBmPMfuuP7QzwLGe3xkmfV0QyCH2Jf2WM+R+rudmf223elnxua75a4F1COobeIpLuMk7jHNb1cwkdOeOe2zHvzdZxzhhjTgFLaI5n9lOQtNYLSCekwMnhrALr8iSO3wPo5Xi/htBZ70nCFXxPWO/zCVc6rY9jzmzClZkxzUXof5AqQsqtTOt9nzjm7e94/8+EztoAlxOuuNtBSJEY1+/CWv8vgYUR7c363D7zNvtzA32B3tb7bsDvgSnAy4QrUb9tvX+AcCXqcr81xTFvf8e/yUKgJNl/Y60uLHz+USYT0qBvBx5O8tiDrV/Qn4Et9viEzp9vAVutn30cv4B/t9byITAyxvl+TWjbXEdIyn8rnrmA+wgp1LYB0+Kc97+scTcBpRFfrIeteT8GJiXyuwCuIbT93QRUWK/Jzf3cPvM2+3MDVwJ/subYDPzQ8fe23lr/y0AXq72r9XmbdX1wtDXFOO/b1jNvBl7grKUmaX9jYt2kKIoSM6mqA1EUpQ2gAkRRlLhRAaIoStyoAFEUJW5aXIDEE5ilKEpq0qICRETSCJmPJhEKGLpHRC7z6T/d61pz01pz6zN3jLnbyzO39A4k1gClVvtHbsW59Zk7xtzt4plbWoC0RJCcoigtRIs6konIHcBEY8z/sT5/A8gzxnzH0Wc6loTs1O2cq9PPDQVrpncScvufQ+W+L6g/YxrbBvXpTo8u6URj857Pw6KCOomQneV974EDB+jbt288j5kQrTVva86tz9y6c2/cuPEzY0xci4n+zUsuUYN1jDHPAM8AdOk/1PQvWNh4rfM5nen7xemwAU8CS2eMIS8ny3fi7MKyJm0nrPvn3ppLwdjBMTyGorQfROSv8d7b0keY94GhVoq3zoQCiEq9OndND1/e3i9O0zlixQa4c/Fa7n/hfd+JJ13Rz7XdAEWllWQXlnHVj95gfdXB6E+hKArQwkcYACsX5UJCEY/PG2Pme/UdOXKkOffuJ9j+2Ymw9t5d0/j8ZEOTRAXndk3n2YKRnruRXYeO851fbaRij19ailCkke5KlI6CiGw0xoyM695UDqbLHTbcXDlzEbsOHmVvxNElq3s6aZ06UXM0vD1dYOgFvSi+9fIwQbJ0zQ7mvVJJg4F/GDeI7KweFJVW+s7fLaMThZMuUUGitGvarQDpMeBi0/cbPw29zxCO1TVd65DzujXZoQB0SYOP5+c3fh5cWMaZiD6TrujHtHGDmfXrP7H3i1Oe66guyfe8pihtnUQESEq7stvWFoBjdYYh53Vr0mf7ZydYPmMMgzLDr51qCO8XKTwAyjfX8NSbH7FmzniqS/IpvjWXdLekboqiuJLSAiSze0bY5x2fnXBN2njn4rV868vZLJ8xhi5pobbcC3qE9bnuYne9yLqqWrILy7jk4dcA2PZ4PtUl+QzKDOWftX8qitKUlD7CjBw50qR99XH2Hzmr5+jXszPpnWiiEwEojqL43HXoOBN/+juO17ntR0JMuqIfi6aOSmzhitKGaLdHGICn7x0RtuuoOXqahfe4F4aLphQd2Kc7b/zzdQzK7OrZp3xzDXctfi+epSpKhyPlBUheThYvzRgT1jbzVx9QfKtbtnwYPvd13/EG9unO6gdv8rwfQseaoXPKWLpmR+wLVpQORMoLEAgJkfN7dW78bJtuq0vySYtQitSebGDonDLuf2EDuw4db2zfdeg49//XBkbPX8Wrm/ZQMHYw1SX5LJ8xhgyXf4W6M6EdzXVPvK3OZYriQZsQIBA6yjgpKq2k+JUP+fX0MU361p2B8s37+fqzf2xse/4PVZRv2c/+I6eYuayicXeRl5PF27Nv8DzW/PXQCb753NokPomitB/ajADJy8li2rhBYW1L1+wkLyfL00195+GTnvqMBW980vjePtaMzunt2vdkfSiWRnUjihJOmxEgAEW3DAsTFgVjQwJl0dRRns5e66pqKXh+Lfddk9O4y0gTmD3x4iZ9X5oxjuqSfN+xil/5MNHHUJR2Q0tH4yaMn4m1+NZcV0vM7z45yPT/2sCCO4dHjdq1GZTZhZ2Hm3qnLnlvJ9lZPdS9XVFoYzuQaBSMHexpXancd4Rpz68PPNbqB8c3cUaziWYuVpSOQrsSIBASIl+72j3J2TEfBzI3ymddz9/0aeo+D6iJV1FIcU/UC4ZcZrre8STn9czgP75+NXk5Wcx++U+s2Lg3rC2SXYeO8/Vn/8jOwyc9x542bhBFtwyLuob1VQe5c7G3FSazezqLv+GdQkBRUp12G43rzEiWe0EvymddG5ZZrGs6fPToWYXn0jU7Go8XtkVlXVWt5/h+UbbOeaaNG8RvNu6l9mS9Z//cC3pQPut6/wdSlBSkXbuy22yrOdKk7WQ9YaZVp2l2XVVto1Ul0vxr45bm0I3/WrOTirkTqS7Jp3Ok55pF5afHonrBKkp7I6UFiPOrWnfGXe/g3GE4TbPOB/M7qgTx76h3bNLGDvE+qtSebCDv0Td9x1KU9kRKC5D+54Z7h857pdLVacx2NbetML26pFHkE+sSybqqWl+l6PnnnHWjn3f7MNeUAjY1R+ui5mdVlPZCSguQrJ5d6Jx29nO9gWnjBjfxGL1r8dowIfJh8c1N/DT8gucgZJp1ChH72HP+OZ152hH9O7BP9ybBfZGUb64h9wflGkOjtHtSWok6cuRI03fqT6j89Kz+w02ZCtC7azoVcyf6jrfr0HG+/MQ7vn2Cpi9cX3WQ6Us3+CpWBajSdIhKitOulajFt10e9tlNmQpQe7I+qm/GwD7dWR5l9xBUsZqXk9WoWO2R4X6oMcCQgOMpSlsk5QVIXk5WWLi9rUx104UE8RDNy8mKusvILowtF8iWeZM9hUgD0XOUKEpbJSEBIiLVIvKhiFSIyAarrY+IrBSRrdbPTKtdRORnIrJNRDaJiHtaMRcemRKuvygqrWTR1FGu0bNBg92iCZGi0sqYom+3zJvsmT+19mQD2YVlqlxV2h3J2IHcYIwZ7jhDFQJvGWOGAm9ZnwEmAUOt13RgUdAJCsYOdrV8vDRjXJOEQkve2xl44X5+HeDvhObG6gfHhyU+iqR8c01M4ylKqtMcR5jbgKXW+6XA7Y72X5oQa4HeItI/6KBzPawoXx3hHvcSlJ/c9aWE7o/k6XtH0LtrmwtyVpS4SMgKIyJVwGFC+sLFxphnRKTWGNPb0eewMSZTRF4FSowxf7Da3wIeNMZsiBhzOqEdCoMGDbr6r389W/fXqeDskiY8ddeX+NJFma6WlU5AkSNL+6ub9jBzWUXjdWcG9yDWGWfsTPErHzbZ6WR1T2fjD89aga798SrXdADDL+rF/8681ncuRWlJWtMKM84YM4LQ8eQBEfH7ZridFZpIL2PMM8aYkcaYkX379g275nRJP9Vg+MFvtjCwT3dXheoZQnoM2xdj3it/CbvuVLgO7NM9qk5kyXs7GxWrv3A5Jh08Hm7OXf3geNdxKnYf0Uhepd2QkAAxxuy1ftYAvwHygP320cT6aR/8dwMDHbdfBOyNZb6iW4a56izsjGRu6gzbyWzE32RGHd8rZqZx/tJKXt20h34+eg4nxbfm0i1DSItoLyqtZPbLFWFJnxWlLRK3ABGRHiLSy34PTAA2A6VAgdWtAPit9b4U+KZljRkDfG6M2RfrvD+560tkdksns3sG874a7iNyUWbT3B0G+OZza5kz+bKwdrdkQUW3DPNNaQgwc1kFT987orECnh8FYwdTOW8y213GW7FxD+Of8j82KUqqk4i273zgNyJij7PMGPO6iLwPLBeRbwE7gTus/q8Bk4FtwHFgWjyTTrlyAFOudFecPnnHl1xzd5yshy8/8Q6jc3rz0oxx8Uwbxp2L15J7QQ8euHFoo16lWxz/kqcaYNLCdzUNgNJmiVuAGGN2AE1MGMaYg8BNLu0GeCDe+YKQl5PF8hljuGvx2qbKFWIzy1aX5DOksIwGj+uVnx7zFWaRZAjUuSyq8tNjgdekKKlGynuixopdyS7IESMa20vyfSNvY2Hr4xoTo7Q/2p0AgZAQ+Xi++xc2u7CM7MKywB6rXv4n8YwVTUnrRfErHzbOlV1YxuyX/xTXOIqSbNqlALHxKhQF0T1Whz4U+rIWlbrnIIkcK0gQXtEtw6KmFfAa38mKjXs1cZGSEqR8OP+GDRuidwyA2xc8Q7yPFs7+3TKEynmTAcgpLHPVr9jMGj+EWeMvTWitXuuIJGjqAUXxo12H8zcnbkpNN07UmcbAuu+OH+Lb9+lV213b7R1NdmEZkxa+G8syPdEoX6W16TACJLO7u8Fp8EPBQvdtC0603UUD7hHBTmEV1PJyycP+x6Lak142IkVpGTqMAFn8Dfcd2hkT7vLuhVOf4hdxCyGdRSLpDJeu2cHFD5dxSuWDkuJ0GAFi+4h44cyrCjR6pNovpwPa0/dGT2Vy5+K1vLppT+Nnp+ern1J2fdVBikorOR1AePgpiRWlJegwSlSbcSWr2FPbNEoWYsthGq1inU2sis6hc8rwqsDZCSgIWFFPUYLSbivTNYcA2XXoODc++Y6nAtXPMuNFNBNuUCHiN45WvlOaC7XCxMDAPt3Z+ni+azAdhJSdfzunLCYdRnVJfpOIWydB0hn6JV9W4aGkKh1OgNiUz7qefj0zXK/Vn4F7nlkbkxDZXpJPuo/fe/nmGk+P1UkL3/WMucnslq7CQ0lZOqwAAVj/yARP9/IGA/c+E13H4eSdf73BN3Ym0qN016HjFJdu8TXrLv5mXDtLRWkROrQAgZB7+deudo+orY9RPWRXrQsagDd7eQVL1lR7Xi++NZe8HO9avIrS2nQ4Jaob0RSrAGkCP7wlN6xk5vqqg3z9mbXUGRiU2SUsjaGfQrQT8LN7h/PdZRWuR5eMTqFSFpHlOddXHWTmrz6g5ujpxrZpapVREkSVqAliK1a9dCIQOtLYKQ1tin67pVHoRCZQ9rO8nCGU2SzTZT4Btj6W30R42PM5hQeEjkWaHlFpLVSAOFj/yIRAKQ0biXJW8bL02Hx2tK5JCQjfWBuP+VZs3MONT2p6RKXl0SOMB15lGWxG5/TmexMuDXMm8xM8fkeafj0zWP/IhKhrWl910DPbWizjKIoTPcI0A6sfHO/rKr6uqraxzm60XUs0ao7WBeqXl5NFVUk+HmV4A4+jKMlCBYgPL80YFzWZUFCiOZvFwtbH8+ndNVmjKUr86BEmAJMWrqby0yNN2oPsOu5/4f3Gmrg9MoRjPqaeWN3oveJmJl3Rj0VTRzVpX7pmB8WllZzx6aN0PJr1CCMiz4tIjYhsdrT1EZGVIrLV+plptYuI/ExEtonIJhEZ4binwOq/VUQK3OZKVYpvu7yJ/rKH1znCQXZhWVhB7WN1xjcvap0J3WMnL4rGI1Pc0yN6FfEusoSH3UfTIiqJEnUHYpWrPEqoMPYVVtsTwCFjTImIFAKZxpgHRWQy8B1C9V9GA/9mjBktIn2ADcBIQrWeNgJXG2MO+82dKjuQeMh79E1XnYRz1+IXeeus3RuEYT98nSOOHABuuyM3Ra7uRJRm3YEYY1YDhyKabwOWWu+XArc72n9pQqwFelvlLScCK40xhyyhsRK4OZ4FpzpL1+wgu7AskEJz/GXne14rKq0MnO0d4JqLzwvc10n55hqt1avETbxK1PPtspTWT1vTOADY5ei322rzam93zC+rjN7JYs7ky3yVodEyx0eONXzAOXTCu3yEW+1gCC80riixkGwrjNufqPFpbzqAyHQR2SAiGw4cOJDUxTU3lzxc5ptJLNKiM7BPdyrm+m/EgpSLsMf63+98mR0l+Z6u7T+8xb/GjSZpVmIlXgGy3zqaYP20tXa7gYGOfhcBe33am2CMecYYM9IYM7Jv375xLq/lyS70zmGa0Smk0/DSNUSz5mQXBkv8HI2CsYOp9vEjqT3ZkLSM8UrHIF4BUgrYlpQC4LeO9m9a1pgxwOfWEecNYIKIZFoWmwlWW7vAb5cwKLOLZ2yLk2hCpKi0MmpSoqD4mYorPz3GtT9elZR5lPZPEDPur4E/ApeIyG4R+RZQAnxFRLYCX7E+A7wG7AC2Ac8C3wYwxhwC5gHvW68fWW1tFltZGi0NoTNCNxo/v3e47/VkKjz9KuTtPHxKFatKINSRLA6WrtkRVfF4btd0/jx3Ysxjv7ppT3jAngs/v3c4U65MXAdd/MqHvoraaWOzue+aHAb26Z7wXErqokmVW5ggis1EvuTNlfE93rnUV6R9o8F0LUiQrX3xrbkJ7RDsIL2WIC8nK9DRSVHcUAESIwve+MTzWreMTiyfMcZVYZr36JuNOpOgytAg1hlngqN4mXLlAKpL8n3Tm6jbu+KGCpAYmT3xYs9rlfMmeeYwdXqmlm+uCZzxPZoQmbmsImkKT7+iWjVH65JmBVLaDypAYqRg7GBfC0ZQ7ly8NvAX36swuE1RaWWLWE3U7V2JRAVIHNgOWdPGZic0TlAX8sXfGOnp/BXrWNGIlv9E3d4VJypAEuC+a3Iav9i9Osf3TxkkfD8vJ4utj+f7Fge3x0pUV7Fo6qhGIeKVsiAZ8yjtA/+9seKLnc09CLkX9PAsILWuqjbQGLZ1xs+MnIy0hk6T7dCHylzLXWj6RAV0B9JilM+6niHndUvKWMlMjxiNrY/7W2eUjo06krUwft6fkcWpouG3Exmd05uXZowLa3N60MaaPtFtrljXq6Qm6kjWhii6ZRi///4NrkWs/MpIuOFnDVpXVUvB8+EepiXlHzW+96vC54ZbjpGdh0+RXVgWU+IjpX2hAqQVGNinO+sfmcDyGWPo4jiLDMrsEtM4tjXIy1fkd58cDPPd6Nera1zrhZDg85pryXs71bzbQVEB0ork5WTx8fyzdWUSOQ547UbKN9c0eqs+eceXGgVWtKp5sVJUWhnYOU5pP6gVpp1QMHYwWT27uEbyzlxWwcGjpygYO5iP5yceY9M1XThZ3/QM9M3n1vLRoy0Tw6OkBroDaUdMuXKA506kqLQycLmIaPzyW6NJd/nLOVmPJiPqYKgAaWcUjB3s6U26rqo2KSkL83KyWPaPY1zNu7ZiVeNmOgYqQNohi6aO4rqL3YP6Kj89lhSrSV5OFi9FKIGdaAqAjoEKkHbK0vvG8LWrL3S9Fku5CD9sJXCs1iOl/aBK1FZkfdVBvvH/1jbJ5u7mBBaNSAe1aeMGseCOq+jVpTNL1lQ36W87hsVaAQ/C3dvtvK9ujmbD574etWyF0rbRHUgr8tSbn7iWglhXVRuzDuEXEbsKW5jcd00OgzK9/T/iyfbudEKz43vc9C61JxtUH9LOUQHSinxvwsW+OoRYdBX9enUO+5xh/WYH9unO6gdv8k1M5PQViZdFU0d5OpolY3wlNVEB0orYOgSvYLUl7+0M/MV7+t4R9OockkZpAo9MaWrO7emTcmDmsopm/ZLPXFahjmbtkCB1YZ4XkRoR2exomysie0SkwnpNdlx7SES2icjHIjLR0X6z1bZNRAqT/yhtF79UgtFKPNjk5WTx4Y9uprokn+2Puxeyen5anm9kbdC57LgYu+JeUIJkmlfaFlGjcUXkWuAo8EtjzBVW21zgqDFmQUTfy4BfA3nAhcAqwE4i+gmhIlS7CRWXuscY8xe/udtjNK4X1/54lW8wXTzKTi+ilaVINCP82MdWsveL067XYo0CVpqfZo3GNcasBoJWkbsNeNEYc8oYU0WoQl2e9dpmjNlhjDkNvGj1VSxWPzie0Tm9Pa8ns7Rlc2d7XzPnK67RuxBSwA55KDnZ5JXWJxEdyEwR2WQdcTKttgHALkef3VabV3sTRGS6iGwQkQ0HDhxIYHltj5dmjPP84kFIGdlSQmTmsoqEHM6KbhnmebxpMImPr6QG8QqQRcAQYDiwD3jKanc7Yhuf9qaNxjxjjBlpjBnZt2/fOJfXdvH74kHs1hk/ls9wd0e3WfLezoTmipbBXtMAtH3iEiDGmP3GmAZjzBlCRbTzrEu7gYGOrhcBe33aFRfsPB9eLHlvJ9mFZVz1ozcSsmzY7uh+2HPFG0NTMHZw1KOZJiVqu8QlQESkv+PjVwHbQlMK3C0iXUQkBxgKrCekNB0qIjki0hm42+qr+ODlim5z+Hg9T73pXSkvCHk5Wfz++zdE7eeVEDoIL80YF9XdPdIRTmkbBDHj/hr4I3CJiOwWkW8BT4jIhyKyCbgB+GcAY8wWYDnwF+B14AFrp1IPzATeACqB5VZfxYcFd1wVNfHP5GHnJzzPwD7dm70W7+oHxzP8ol6e19UhqW2iSZXbAKPnr2T/EXezKIT8MopuGZaUuaKZeCdd0S+s7EOsRHuWRMdXYkeTKrdzVtw/zjeepSWVkYm6pT9974hmHV9pWVSAtAHseBa/ynRFpZVc98TbCbuLB/EsTcQEm5eTFbXOTDILhivNiwqQNoRdmc6Lvx46wfSliR35nBnNvEpbQmI5RebemkuvLmlRrTNK6qMCpA3iZ9GoPVmf8Ph2ZO2WeZN9+8WbY7Vg7GA+LL45qnVGjzKpjwqQNsjqB8d75j0FuOKHrzP+qXebHGfWVx0k79GVZBeWMXROWaBjgt+OZ11VLeNK3moyz6ub9jB6/kruf2EDuw4dj/osXt63s16siHq/0rqoAGmjLJo6ylOIHD3dwLYDx/jmc+HRr99Z9gE1R0MWkLozoWNCokJkT+3JJlG2M5dVsP/Iaco37+fGJ9+JOn7RLcM4PyKfCUD9GQLdr7QeKkDaMPZRw0vxebI+lFbQpsbFfLrgjWCOaH47HvA2/wYtofn0vSPo3bVphs06AxfP8TctK62HCpB2QMHYwZ5xLbUnG8ixomv/weWoMHvixS53NcUv45hNIlG8eTlZVMyd6Dr+6TOQ9+ibcY2rNC+aVLmdYMe1uCVpNsA/LavgxRljkuJwlga4pHIFQseXLmk0rsHP0uI5voQidp3UHK3j/hfeVyezFEN3IO0IO0Wim/v7GUIZwZLhX7G9JJ80H2eOUw0hf5LqkvyYs8sDfHnoea7t5ZtryP1BuaZGTCFUgLRDymdd73ktWf4V/3bPcN/ricwz7/Zhnp63J+rOcJemRkwZVIC0U/wSEw2JEu8ShClXDvD1jIXocTVe2J63Xo5sqRu91fHQYLp2zNI1O6LuBM4/pzNP3zOCvBz3UphuLFz1Ef+2ajsGGNC7K3tqT/r2jzXYb33VQWb+6gMOHD1NJ9z1Lbn9e1F86+UxrVtxR4PpFFds64xX7RmA/V+cjtn9faElPCDkBxIttiVWt/eZvwr5qxhCwsPtj7Ry3xE9yqQAKkDaObZidfmMMaR7/LZrT9bH7ZYOoSNFNGtLLMcZ29nNpsjDz8VY4ybjSKbEhwqQDkJeThbL/tE7B+q6qlqGPhTsi9gtwvjfu2saL80YF8hPJIgVKCPir9JO8eg1dgPhDnNKy6E6kA7IkMIyTz+OC8/pzJo5X0lo/ESTEkXqbpyCw2/sIed1463ZN8awUgVUB6LEyPaSfHp3dVeM7P3idMKlI6JZZ8o31/j6ctjZ3Ht1SWvipu8WM2Oz/bMTScmJogRHdyAdmJsWvM32z064Xsvt35NnvjGKgX26xzX2q5v2RC2V+fN7hzPlStfyQJ6srzrI9KUboqYt+P33b4h77R2NRHYgKkA6ONc98TZ/PeQuRAZldmX1gzfFPfauQ8e58cl3fAPqEknmfNfi91hXVet6TUtoBkcFiBI366sORi16nXtBD1/v1mgEscB87eoLWXDHVTGPPfShMl8BlQ5sa+aM822dZtWBiMhAEXlHRCpFZIuIfNdq7yMiK0Vkq/Uz02oXEfmZiGyzSl+OcIxVYPUOG9M1AAAWiUlEQVTfKiIF8SxYSS52XRif7IVUfnosIb1IkF3Gio1744rT2fp4Phee460XqQcuefg1zbHaTARRotYD3zPG5AJjgAdE5DKgEHjLGDMUeMv6DDCJUEGpocB0QmUwEZE+QBEwmlAluyJHTV2lFRnYpztbH8/3FSLlm2sSmqO6JJ/0KB5nQXOTRLJmzld885WcajBxj634E1WAGGP2GWM+sN4fIVQYagBwG7DU6rYUuN16fxvwSxNiLdDbqmQ3EVhpjDlkjDkMrARuTurTKAkRTWcQNA2iF+/8q38FvCOnGuLe6SyaOopePi63R041qHWmGYjJjCsi2cBVwDrgfGPMPggJGcD+L2AAsMtx226rzatdSSH8Ngl1Z2DeK/FH2Q7s0903yA8S2+k8/vf+8TZ3Ll6bsIlaCSewABGRnsB/A7OMMV/4dXVpMz7tkfNMF5ENIrLhwIEDQZenJIm5UerC1Jv4s7EDgYLq4o3iDWISTvQopoQTSICISAYh4fErY8z/WM377SLb1k/7N7MbGOi4/SJgr097GMaYZ4wxI40xI/v27RvLsyhJwOnE5RVOH4vbuxvXXRw9gja7sIxJC9+Neexp4wbRSfxr2lzysMbOJIuoZlwREUI6jkPGmFmO9ieBg8aYEhEpBPoYY74vIvmECmlPJqQw/ZkxJs9Som4EbKvMB8DVxphDXnOrGbf18XN7Xz5jTMLh9H7jg7sFZ+maHcwvq0QQ5uRfSsHYwa733v/C+547juYuJt6WaG5X9nHAN4AbRaTCek0GSoCviMhW4CvWZ4DXgB3ANuBZ4NsAlqCYB7xvvX7kJzyU1MDP7T3RKnj2+H7pEd2YX1bJ6YaQdcUv38miqaNiHluJDXUkUwIRTS/h9PwcPvd1ak+e3VdEC9AL4vZu7xgmLXyXyk+PuV5zwyupUqLOce0JDaZTmp1odWHqDFxV/Aa7Dh0PEx4QCtDzY8qVA6IW9c4uLKP4lQ+bCI9+PTN87ysYO9i1T+Wnx7j956t971WiozsQJWb8gvDSO4UqykUSi84hFitMouMW35rrqUPpKOgORGlR3pp9o6c/h5vwCGJ1cTI6p+UclItKKyku3aI1eONEBYgSF0W3DGNA7y6e1wUas4gtvc8/P0gkC+7wLxlh41b/Jh6WrKlm/FNagzceVIAocfPi9LGe8TMGuPbHq+Iad2Cf7oGOJrEqQf28YCOr+SnBUAGixI0dhDco030nsvPwqYTGry7JxyehfMwU3TKM6hL/oEElNlSJqiSFvEffpOZonef1WBL8BKlnE0mQ2jOzX/4TKzY2cX5uxLnrcfbt1zOD9Y9MiGk9bQlVoiqtzvpHJviaev2S/jjJe/TNuMpi+tWeWbpmB8OKXvcVHgAFz59NrOTsW3O0LuyachYVIErSWDR1FNUl+Qw5r5vr9b+dUxY1pN5tF+NUyPrhFWk7v6ySIy5Kjsgxf/fJ2bV1iShd8btPDnLTgrd95++IqABRks5bs2+kk4ueof4M3PPMWl8h4nT6yr2gB9Ul+VQF9PUo31zDq5v2uFyJXelxyiVn8/bPTlD8yocxj9WeUQGiNAsFY90tHg0G7l681tPvYv0jExp3Bm5WlmiiYOayiiZC5OH8S+nVJY3MyIpYFrPGD6FzmjBr/JAoo4eOSrNfrlC/EQtVoirNxuyXK1ix0W1HEKpmVzE39oR0QRWsbsedXYeO8/Vn/8jOwye57uIsX/+UWGJ/2jqalV1JSaKVdUj0SxjtS55IyH7xKx9GLQreXqwzKkCUlCaaiRfij0kJEjeTaLxLtNIRbV2QqBlXSWnWPzLBM6eITVFpZVxJj4PsMopKKxPKhbr18XzSfb4pNUfrOmzpCBUgSotQMfdmhl/Uy7fPU2/GV3ph+YwxUZWrieZCXfaPY1wtSzYdtXSEChClxfjfmdeSe4G3EFlfdSiu/8XzcrJ4KUpBbyChqNu8nCxenO4/x7FTDR3OOqMCRGlRim+73POagbi8UCH0Bf/5vf5RvEvWVHP9E/FH3Uab4wwkNH5bRAWI0qLk5WRFPW5kF5aRXVjm6z5uu6c7dyxByjo0AMN+WB7Wtr7qIHct/qOrDmbSwnfDMsRHm6OjBfWqAFFanKfvHc75vTpHTUfodC2PZH7ZRxw51cD8so/C2qMVrgI4cjo861HRb7ewruoQRb/d0qSvnULRmUoxSOmIjoK7a56iNCNTrhzQ+D+5nxnWq1Ll+qqDnG4I2VXtnzZFtwxrjMoNmhrxwJGTYT+j4Zxj6Jwy6lyysHUUVIAorcqkK/p5WkhONYQLgdE5vfniRF2TxMp/O6eMZf/YtEaN39j2uNUl+Xx2LOSjYv+EYI5kAI9MyW2it8kuLOswWd+jHmFEZKCIvCMilSKyRUS+a7XPFZE9EbVi7HseEpFtIvKxiEx0tN9stW2zilEpHZxFU0cFOnZAqCJepPCAUJCe2/Fj0dRRUfOxOgWU88sQRHhAKOv7eT3cs753hDq8QXQg9cD3jDG5wBjgARG5zLr2U2PMcOv1GoB17W7gcuBm4D9EJE1E0oB/ByYBlwH3OMZROjB2pjD7tTyASbYJHuqIpfeNCSSgOgH/5AimixxudE5vz3t/88A1XDf0vCbt5Ztr2n30blQBYozZZ4z5wHp/BKgE/FTRtwEvGmNOGWOqCFWoy7Ne24wxO4wxp4EXrb6KEkZeThbVJfl0jvLXaeswM7unU3yrt3nYFlA9fQY8A/z9iLOCJi8iM/xLM8Z53juwT3eWfmu067Ul7+1s174hMVlhRCQbuApYZzXNFJFNIvK8iNj/4gOAXY7bdlttXu2Rc0wXkQ0isuHAgQOxLE9pZ3zymL+b+tbHQzuWP/1wYqAavc9Py/M1Id/803cb3y+4YziDMrsCwctSeH2ZbmjHviGBBYiI9AT+G5hljPkCWAQMAYYD+4Cn7K4utxuf9vAGY54xxow0xozs27dv0OUp7RS/o8Po+as8Egi5U7nvc3p6mXaAY3WmMdfHwD7dWf3gTRTfmssHf60N5CFb5FFdrx7a7S4kkAARkQxCwuNXxpj/ATDG7DfGNBhjzhAqop1ndd8NDHTcfhGw16ddUTx5acY4Tx3G/iOn+OcX/Wvqwlmns7mlodSG6T7bkBUb9zDhJ+82fuGLrHuCeMgWjB3M3/RxT+d4w5PvxBUsmOoEscII8BxQaYz5iaO9v6PbV4HN1vtS4G4R6SIiOcBQYD3wPjBURHJEpDMhRWtpch5Dac/YOgw3QRLEB2PBG59w5FRD43a33vgXpTpRb/j6s3+Ma61P3vElendt6h1Rb+Cbz7W/xMxR84GIyDXA74EPCemaAOYA9xA6vhigGphhjNln3fMwcB+h3dssY0y51T4ZWAikAc8bY+b7za35QBQ3YsksBqEdyII3PuGyC3vxfnUtBWPDS0B4OZzlXtAjzGwca4Kiy39QzjEXCTfpin4smjoqprGaE00opCgJ4Oc01q9nBjVH6xid09vXEuPG+KfeZduBpn4rkFi2tGSjCYUUJQGKbhlGsYcCtOZoHdUl+TELD4DH/m6Ypzt+e0F3IEqb4v4X3vd0T4+2S4jcaUQmdi4u3cKSNdW+8/u5qE9a+G7jkceZ79XtiBRvUunmIJEdiMbCKG0Kv8xi66pqXdudX2wntScbfK+74dbPLVN8nYGFqz5i4artruPUnmzgpgVv89bsGwPNm6roEUZpU/iVz+wErqbSaMIhqPCwifQ98TLxPh0hPCL1Hts/O8FVP3ozJl+WVEMFiNKmWDR1FMtnjHHVLZwB7ly8NvAX0stJLaMTvm7vbsWr3OgeQAFy+HgdM5e13UJVKkCUNkdeThYfz/e2YsxcFu5cFunzMSizS5hi1Olf0jlNeGRKLsdO+zuYOOfw8il5/O/PmoqjfdFufLJtururElVpswyf+3qjHiOSRItWBc0H4jyWOO/xUpL6JTlqrWp36geidFiCKEATLSyVzOJV0YpUOZ3XBmV2YfWD4wOvM17UD0TpsJTPup4h57nHn9gUlVYmVPQpaPGqIMWl7AhiL2WwUxjuPHwqtoW2AipAlDbPW7Nv9HQEs0m06FOQ4lWxFJdaNHUUPbu0/a9f238CRSEUCTttbLbn9SOnGhIuLBWkeFUsxaWe/4c818C7toQKEKXdcN81Ob7Xl6yp5stPvMPsl/8UeMziVz4ku7CM0Y+tBIhavOoMcNOCsxYVt/o1dy1+j+zCMu5cvJZL+vcMvJZURAWI0m4Y2Kd71FozACs2Bk9DY1tV9n9xmqLfbglUvMppAXarX+P0mF1XVetpBvZLOZAqqABR2hU5ff0LeCfCp5+fAKBbRrCvzdI1Ozzr1zjxiq1pC2Uh2vYBTFEi+N6Eiyn67Rb2fX6C2hP1nv0i681EBuG5Be0dPlHPXYvfo3DSJcwvq+S0Tx1LL9Nv0GJXzr7TxoXnL0kldAeitCvycrIon3UtFUUTqS7JZ1Bml6j3rKuq5fIfvBbW5hW0t66qloKxg/lkfsgcG836Y+N1tLLbv3b1hZ73LnlvZ2Nt3lRDdyBKu8bpiOUWNWtzrM7w6qY9jToOr6p2kfEzBWMHNzqQ+e0u1j8yoUlbmqN9wR1XseCOqxqvRY4Va8BfS6ECROkw2F/2i+eU4Rbq8k/LKhoFSDwpB9MEfFQdTbjIIwFzW0Jd2ZUOSTRdhJ231NnPSxexdM0OflRaiY9KJGxc587GKwmS1/riSa0YDXVlV5QYiWYiLd9c08RfxCu4bn5ZMOFhj+vEMwmSh6v7uqralCqXqQJE6ZCUz7q+sRavV+W5SH8Rb+ttuJN7t4xOCftw2HlP3NKSBC383RIEqQvTVUTWi8ifRWSLiBRb7Tkisk5EtorIS1atF6x6MC+JyDbrerZjrIes9o9FZGJzPZSixMLS+8YEsqY8MsW9z8P5l9IlHbqlC5ndM3jyjispn3U9oyPq68ZKXk5W1PKerU2QujAC9DDGHLUq1P0B+C7wL8D/GGNeFJH/BP5sjFkkIt8GrjTG/F8RuRv4qjHmLhG5DPg1oQp2FwKrgIuNMZ67P9WBKC3J7f/+eyp2fdGkvXfXtMa8I726pPHcP4xyrcUbmUPka1dfGMjr1a1OjJ2mQHCp/2rd84etB5k98eKEUhVAM+tATIij1scM62WAG4EVVvtS4Hbr/W3WZ6zrN1lC6DbgRWPMKWNMFbCNs+UwFaXVefqeqxsLajtxJi06cqqB6Uub/qc2aeG7TY4WKzbuJff86LEutl7k/hfeJ7uwjPtfeL/RbOv133v55prAJTebk6C1cdNEpAKoAVYC24FaY4zt6rcbsIMEBgC7AKzrnwNZznaXexSl1bELaleX5JN7gbdLfO3JcA/XoQ+VefppVO4/6toeidPz1S/zfKoRSIBYRbSHEyqInQe4HQZtYemWNsH4tIchItNFZIOIbDhw4ECQ5SlK0im+7XLf/B9DHzprZvXLMBaURITG/S+83/h++NzXyS4sY/jc1xNfVABissIYY2qBd4ExQG8RsR3RLgLsw95uYCCAdf1c4JCz3eUe5xzPGGNGGmNG9u3bN5blKUrSyMvJoqok39MFvc6cFSIZEZKmd9c0qkvym4T+R/ZLFuWbaxpNu/ZxyytXbLIJokTtC9QZY2pFpBvwJvBjoAD4b4cSdZMx5j9E5AFgmEOJ+nfGmDtF5HJgGWeVqG8BQ1WJqqQ6fsmbAdIFlk0f46pYjSSIA1vkbiS9E2x7zLvKnU11SX7Y9R4ZwpZ5k6OuqbkdyfoD74jIJuB9YKUx5lXgQeBfRGQbIR3Hc1b/54Asq/1fgEIAY8wWYDnwF+B14AE/4aEoqULF3Jt9867Wm1A9mpzCMu5avMY3I9ms8UPwqxbjdpRZePfZncz5vToHWjOE4nsWrvooescEUFd2RQmIXzCek3Rgm08i5rGPrWTvF6cDz1t8a27jvP16ZnC63jRR5ELTHYiz3Q91ZVeUFqBg7GB+//0bovbzzkISIhbhMfyiXswvOyu0ao7WuQoPgNkv/6nZ9CxeqABRlBjYZ2Ulaym27DmCuwGzKSs27uXyAc2Xkc0NFSCKEgNPvRmsbEN2YRm5P3iN+1/YEKYTiQzQ6+FTgxdC1h7BIITEiF9xcYCK3UdiGj9RVIAoSgx8b8LFgfueqDOUb97P15/9Y2NbpGt7tBq8AKcaQg5Tf3/1haz6S2z+IqcCjJ8IKkAUJQbycrKoLskPK8id3ilUeCqzm7vPyM7DJwON3Ql8UzCu2LiXuhjlQT0kVJUvGpqRTFHioOiWYU2SC8376uXMXFYR81huhbhjSb4caX2J/FxUWplwwJ0XKkAUJUlMuXIAB4+ecjX1+gmEirk3h5mIY7WktJTbuht6hFGUJFIwdjDVJfkxCwFnTd06Q9gRKRqRXrJulfcis84nCxUgitIMbH08n+UBaunazJ54VjmbIaEjkp0xLVZWbNzb5L5jyYj4c0E9URWlmfHKAu9k2rhB/GbjnsbdRL+eGdQcrUvqOiZd0Y83NtdwhvAkRuqJqigpzCeP5UfNkbrkvZ1hR5FkCw8IxdmccbxPBipAFKUFsJM4pxLrqw4mPIYKEEVpQYKU2mwp7n1mbcJjqABRlBZk9YPjo7qjRxJr/6DUJ0H9qQJEUVqYRVNHxVSY29ZXTBs3yLOGTWuhjmSK0ko4C3MDXPvjVew8fMqzfyoVlLJRAaIorYBbcqIeGdKoaF246iMWrtreGkuLCRUgitIKOD1PbY7VGXIKy8JKFeRe0IO78wY1W/2XSQvfTeh+1YEoSivg9Dx1EqnXrPz0WKN7fGQUcDLwqmcTFN2BKEorUDB2MLn9z+WuxWs9q8/ZOAPxki1AEkV3IIrSSkSrPeNGqilSdQeiKK3M+kcmNGmLLNSdqkTdgYhIVxFZLyJ/FpEtIlJstf9CRKpEpMJ6DbfaRUR+JiLbRGSTiIxwjFUgIlutV0HzPZaitG2c0bhB/UXiJf3c8/8m7nsD9DkF3GiMOSoiGcAfRKTcuvavxpgVEf0nAUOt12hgETBaRPoARcBIQrqijSJSaow5HO/iFaUj4PQXGVfyFntqg6VIDEqnbr3Oi/veaB1MCLvEeIb18tP73Ab80rpvLaEauv2BiYSq2h2yhMZK4GafcRRFiWDC5c3j1h4vgZSoIpImIhVADSEhsM66NN86pvxUROwooQHALsftu602r/bIuaaLyAYR2XDgwIEYH0dR2jdvbklOGH6yCCRAjDENxpjhwEVAnohcATwEXAqMAvoQqpUL7lVwjE975FzPGGNGGmNG9u3bN8jyFKXD8NDkS1t7CWHEZMY1xtQC7wI3G2P2WceUU8ASIM/qthsY6LjtImCvT7uiKAGZcuUAZo0f0trLaCSqElVE+gJ1xphaEekGjAd+LCL9jTH7RESA24HN1i2lwEwReZGQEvVzq98bwGMikmn1m0BoF6MoSgzMGn8ps8Y33YkELf6dTIJYYfoDS0UkjdCOZbkx5lURedsSLgJUAP/X6v8aMBnYBhwHpgEYYw6JyDzgfavfj4wxh/wm3rhx41ER+TjWh0oS5wGfdaB5W3NufeZmQDp375ne+4K/lU6d0vz61X8ev14lpZMqi8iGeJO9ttW59Zk7xtzt5ZnVlV1RlLhRAaIoStykugB5pgPOrc/cMeZuF8+c0joQRVFSm1TfgSiKksKoAFEUJW5UgCiKEjcqQBRFiRsVIIqixM3/B/S67mj32JXHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1472.080010077052"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.spy(K,markersize=1)\n",
    "plt.show()\n",
    "\n",
    "spla.norm(K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 矩阵乘法\n",
    "Kb = K*b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.        ,  -1.        ],\n",
       "       [ 13.48478231, -13.48478231],\n",
       "       [ 52.3417648 , -52.3417648 ],\n",
       "       [  1.        ,  -1.        ],\n",
       "       [ 27.20825099, -27.20825099]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Kb[190:195,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        , -1.        ],\n",
       "       [ 0.07122914, -0.07122914],\n",
       "       [ 0.04526046, -0.04526046],\n",
       "       [ 1.        , -1.        ],\n",
       "       [ 0.04214178, -0.04214178]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 线性方程求解\n",
    "x = spla.spsolve(K,b)\n",
    "x[500:505,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        , -1.        ],\n",
       "       [ 0.07122914, -0.07122914],\n",
       "       [ 0.04526046, -0.04526046],\n",
       "       [ 1.        , -1.        ],\n",
       "       [ 0.04214178, -0.04214178]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LU = spla.splu(K)\n",
    "\n",
    "x = LU.solve(b)\n",
    "\n",
    "x[500:505,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 标准特征值问题的解\n",
    "\n",
    "val,vec = spla.eigsh(K,k=10,sigma=-1)\n",
    "\n",
    "val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sturm数\n",
    "sigma = 0.99\n",
    "\n",
    "LU = spla.splu(K - sigma*sp.eye(nR))\n",
    "\n",
    "dU = LU.U.diagonal()\n",
    "\n",
    "np.sum(dU < 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1651"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma = 1.01\n",
    "\n",
    "LU = spla.splu(K - sigma*sp.eye(nR))\n",
    "\n",
    "dU = LU.U.diagonal()\n",
    "\n",
    "np.sum(dU < 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
